Model comparisons

Libraries

library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages -----------------------------------------------------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(lubridate)

Attaching package: 㤼㸱lubridate㤼㸲

The following object is masked from 㤼㸱package:base㤼㸲:

    date

library(stringr)
library(e1071)
library(ReporteRs)

Loading required package: ReporteRsjars

library(slam)
library(tm)

Loading required package: NLP

Attaching package: 㤼㸱NLP㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    annotate

library(topicmodels)
library(stm)

stm v1.2.2 (2017-03-28) successfully loaded. See ?stm for help.

Constants

theme_sparse <- theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
  panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank())

Load models and data

# Data
load(file = "data/sparse_matrix_no_tex.zip")
sparse <- readCorpus(sparse, type = "slam")
# Metadata
load("data/posts_metadata.zip")
metadata$HW_Tag <- factor(metadata$HW_Tag, labels = c("No HW", "HW"))
# CTM
load("results/lda_final_40.zip")
load("results/ctm_final_110.zip")
load("results/stm_final_90.zip")

Least popular homework topics

metadata %>% filter(HW_Tag == "HW") %>% arrange(Score) %>% select(Id, Score, Title) %>% head()

Shared topic labels

# LDA: data = 10, model selection = 34
lda_gamma <- models40f@gamma %>% data.frame(i = 1:nrow(metadata))
lda_gamma <- data.frame(
  lda_gamma %>% arrange(desc(X10)) %>% select(T10 = i) %>% head(5),
  lda_gamma %>% arrange(desc(X34)) %>% select(T34 = i) %>% head(5)
)
print("LDA")

[1] "LDA"

print("Topic 10:")

[1] "Topic 10:"

apply(metadata[lda_gamma$T10, c("Id", "Title")], MARGIN = 1, FUN = paste, collapse = " ")

[1] "244815 Scaling large range values"                                     
[2] "143414 Combining Standard Deviation (combined measurement uncertainty)"
[3] "175945 PCA on bimodal data: how to standardize the data?"              
[4] "249959 Normalization of data with lots of zeros"                       
[5] "155407 How to compare 2 sets of p-values?"                             

print("Topic 34:")

[1] "Topic 34:"

apply(metadata[lda_gamma$T34, c("Id", "Title")], MARGIN = 1, FUN = paste, collapse = " ")

[1] "198365 Higher SIC and lower S.E. of residuals"                             
[2] "221016 AIC (QIC) - comparision of models?"                                 
[3] "206808 Model comparison with chi-square difference test and BIC"           
[4] " 94292 Comparing nested, non-linear models"                                
[5] " 83042 How can I evaluate my nonlinear model (without a model comparison)?"

# CTM: data = 20, model selection = 46
print("CTM")

[1] "CTM"

findThoughts(ctm_final_110, texts = trimws(paste(metadata$Id, metadata$Title)), topics = c(20, 46), n = 5)

 Topic 20: 
     79511 How to combine two data sets (level-1 and level-2 data sets) for multilevel analysis in R
    58333 What method can be used to test if three or more categorical sample data sets are from the same distribution?
    123621 Can I use K-S test to distinguish distributions?
    52010 Choosing a better data-set
    209410 How can I create a synthetic data set resembling an existing data set? 
 Topic 46: 
     94292 Comparing nested, non-linear models
    194601 Comparison between different models: hierarchical models and model with dummy
    76893 Comparing linear models
    174018 Extract goodness of fit from linear mixed effects models
    215546 Using DIC for model selection: (in)valid comparisons

# STM: data = 20, model selection = 46
print("")

[1] ""

print("STM")

[1] "STM"

findThoughts(stm_final_90, texts = trimws(paste(metadata$Id, metadata$Title)), topics = c(20, 46), n = 5)

 Topic 20: 
     104392 Multiple Imputation: Before or after case exclusion?
    181399 How to combine multiply imputed datasets created with MICE from different cohorts?
    104904 How to replace the missing data from AMELIA results
    228463 perform quality check for imputed data with MICE in R
    230369 imputation of missing values 
 Topic 46: 
     194601 Comparison between different models: hierarchical models and model with dummy
    127353 What model should I use for this?
    221016 AIC (QIC) - comparision of models?
    76893 Comparing linear models
    97901 R fit restricted AR(p) model

Most upvoted post

Topic proportions

# Label the top top_n topics
top_n <- 3
# LDA
lda_popular <- data.frame(
  topic = 1:40,
  theta_d = models40f@gamma[pop_idx, ]
)
lda_popular$top <- lda_popular$topic %in% (lda_popular %>% arrange(theta_d) %>% select(topic) %>% tail(top_n) %>% unlist())
lda_popular <- lda_popular %>% mutate(
  terms = ifelse(!top,
    "",
    apply(topic %>% rbind(get_terms(models40f, k = 3)), MARGIN = 2, FUN = paste, collapse = "\n")
  )
)
# CTM
ctm_popular <- data.frame(
  topic = 1:110,
  theta_d = ctm_final_110$theta[pop_idx, ]
)
ctm_popular$top <- ctm_popular$topic %in% (ctm_popular %>% arrange(theta_d) %>% select(topic) %>% tail(top_n) %>% unlist())
ctm_popular <- ctm_popular %>% mutate(
  terms = ifelse(!top,
    "",
    apply(topic %>% cbind(sageLabels(ctm_final_110, n = 3)$marginal$frex), MARGIN = 1, FUN = paste, collapse = "\n")
  )
)
# STM
stm_popular <- data.frame(
  topic = 1:90,
  theta_d = stm_final_90$theta[pop_idx, ]
)
stm_popular$top <- stm_popular$topic %in% (stm_popular %>% arrange(theta_d) %>% select(topic) %>% tail(top_n) %>% unlist())
stm_popular <- stm_popular %>% mutate(
  terms = ifelse(!top,
    "",
    apply(topic %>% cbind(sageLabels(stm_final_90, n = 3)$marginal$frex), MARGIN = 1, FUN = paste, collapse = "\n")
  )
)

# View top topics
lda_popular %>% filter(top) %>% arrange(desc(theta_d))

ctm_popular %>% filter(top) %>% arrange(desc(theta_d))

stm_popular %>% filter(top) %>% arrange(desc(theta_d))

# LDA results
lda1 <- lda_popular %>%
  ggplot(aes(x = topic, y = theta_d)) +
  geom_bar(stat = "identity", fill = "lightgray", width = .75) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent, limits = c(0, .7), breaks = seq(0, 1, by = .1)) +
  theme_sparse +
  labs(x = "LDA topic number", y = "LDA topic proportions \\theta_d") +
  geom_text(aes(label = terms, y = theta_d + .12), size = 3)
plot(lda1)

# CTM results
ctm1 <- ctm_popular %>%
  ggplot(aes(x = topic, y = theta_d)) +
  geom_bar(stat = "identity", fill = "lightgray", width = .75) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent, limits = c(0, .6), breaks = seq(0, 1, by = .1)) +
  theme_sparse +
  labs(x = "CTM topic number", y = "CTM topic proportions \\theta_d") +
  geom_text(aes(label = terms, y = theta_d + .1), size = 3)
plot(ctm1)

# STM results
stm1 <- stm_popular %>%
  ggplot(aes(x = topic, y = theta_d)) +
  geom_bar(stat = "identity", fill = "lightgray", width = .75) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent, limits = c(0, .6), breaks = seq(0, 1, by = .1)) +
  theme_sparse +
  labs(x = "STM topic number", y = "STM topic proportions \\theta_d") +
  geom_text(aes(label = terms, y = theta_d + .1), size = 3)
plot(stm1)

# Save results
mydoc <- docx()
mydoc <- addPlot(mydoc, function() print(lda1), width = 7, height = 2.5, vector.graphic = TRUE)
mydoc <- addPlot(mydoc, function() print(ctm1), width = 7, height = 2.5, vector.graphic = TRUE)
mydoc <- addPlot(mydoc, function() print(stm1), width = 7, height = 2.5, vector.graphic = TRUE)
writeDoc(mydoc, file = "results/fig/10-comparison-popular.docx")

Similar documents

# Hellinger distance function
hellinger <- function(target, doc_thetas) {
  n <- nrow(doc_thetas)
  target <- matrix(rep(target, each = n), nrow = n)
  2 - 2 * rowSums(sqrt(target) * sqrt(doc_thetas))
}

LDA

lda_dist <- data.frame(
  rown = 1:nrow(models40f@gamma),
  hellinger = hellinger(models40f@gamma[pop_idx, ], models40f@gamma),
  Id = metadata$Id,
  Title = metadata$Title
) %>% arrange(hellinger)
{lda_dist %>% head()}[-1, ]  %>% select(hellinger, Id, Title)

CTM

ctm_dist <- data.frame(
  rown = 1:nrow(ctm_final_110$theta),
  hellinger = hellinger(ctm_final_110$theta[pop_idx, ], ctm_final_110$theta),
  Id = metadata$Id,
  Title = metadata$Title
) %>% arrange (hellinger)
{ctm_dist %>% head()}[-1, ]  %>% select(hellinger, Id, Title)

STM

stm_dist <- data.frame(
  rown = 1:nrow(stm_final_90$theta),
  hellinger = hellinger(stm_final_90$theta[pop_idx, ], stm_final_90$theta),
  Id = metadata$Id,
  Title = metadata$Title
) %>% arrange (hellinger)
{stm_dist %>% head()}[-1, ]  %>% select(hellinger, Id, Title)

How long until one of the CTM/STM posts shows up in LDA?